home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / b1nuG.c < prev    next >
C/C++ Source or Header  |  1988-11-24  |  4KB  |  193 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: b1nuG.c,v 1.4 85/08/22 16:50:59 timo Exp $
  5. */
  6.  
  7. #include "b.h"
  8. #include "b0con.h"
  9. #include "b1obj.h"
  10. #include "b1num.h"
  11.  
  12.  
  13. /*
  14.  * Routines for greatest common divisor calculation
  15.  * Cf. Knuth II Sec. 4.5.2. Algorithm B
  16.  * "Binary gcd algorithm"
  17.  * The labels correspond with those in the book
  18.  *
  19.  * Assumptions about built-in arithmetic:
  20.  * x>>1 == x/2  (if x >= 0)
  21.  * 1<<k == 2**k (if it fits in a word)
  22.  */
  23.  
  24. /* Single-precision gcd for integers > 0 */
  25.  
  26. Hidden digit dig_gcd(u, v) register digit u, v; {
  27.     register digit t;
  28.     register int k = 0;
  29.  
  30.     if (u <= 0 || v <= 0) syserr(MESS(900, "dig_gcd of number(s) <= 0"));
  31.  
  32.  /*B1*/    while (Even(u) && Even(v)) ++k, u >>= 1, v >>= 1;
  33.  
  34.  /*B2*/    if (Even(u)) t = u;
  35.     else {    t = -v;
  36.         goto B4;
  37.     }
  38.  
  39.     do {
  40.  /*B3*/        do {
  41.             t /= 2;
  42.  B4:            ;
  43.         } while (Even(t));
  44.  
  45.  /*B5*/        if (t > 0) u = t;
  46.         else v = -t;
  47.  
  48.  /*B6*/        t = u-v;
  49.     } while (t);
  50.  
  51.     return u * (1<<k);
  52. }
  53.  
  54.  
  55. Visible integer int_half(v) integer v; {
  56.     register int i;
  57.     register long carry;
  58.  
  59.     if (IsSmallInt(v))     return (integer) MkSmallInt(SmallIntVal(v) / 2);
  60.  
  61.     if (Msd(v) < 0) {
  62.         i = Length(v)-2;
  63.         if (i < 0) {
  64.             release((value) v);
  65.             return int_0;
  66.         }
  67.         carry = BASE;
  68.     }
  69.     else {
  70.         carry = 0;
  71.         i = Length(v)-1;
  72.     }
  73.  
  74.     if (Refcnt(v) > 1) uniql((value *) &v);
  75.  
  76.     for (; i >= 0; --i) {
  77.         carry += Digit(v,i);
  78.         Digit(v,i) = carry/2;
  79.         carry = carry&1 ? BASE : 0;
  80.     }
  81.  
  82.     return int_canon(v);
  83. }
  84.  
  85.  
  86. Hidden integer twice(v) integer v; {
  87.     digit carry = 0;
  88.     int i;
  89.  
  90.     if (IsSmallInt(v)) return mk_int(2.0 * SmallIntVal(v));
  91.  
  92.     if (Refcnt(v) > 1) uniql((value *) &v);
  93.  
  94.     for (i = 0; i < Length(v); ++i) {
  95.         carry += Digit(v,i) * 2;
  96.         if (carry >= BASE)
  97.             Digit(v,i) = carry-BASE, carry = 1;
  98.         else
  99.             Digit(v,i) = carry, carry = 0;
  100.     }
  101.  
  102.     if (carry) {    /* Enlarge the number */
  103.         v = (integer) regrab_num((value) v, Length(v)+1);
  104.         Digit(v, Length(v)-1) = carry;
  105.     }
  106.  
  107.     return v;
  108. }
  109.  
  110.  
  111. Hidden bool even(u) integer u; {
  112.     if (IsSmallInt(u))
  113.         return Even(SmallIntVal(u));
  114.     return Even(Lsd(u));
  115. }
  116.  
  117.  
  118. /* Multi-precision gcd of integers > 0 */
  119.  
  120. Visible integer int_gcd(u1, v1) integer u1, v1; {
  121.     struct integer uu1, vv1;
  122.  
  123.     if (Msd(u1) <= 0 || Msd(v1) <= 0)
  124.         syserr(MESS(901, "gcd of number(s) <= 0"));
  125.  
  126.     if (u1==int_1 || v1==int_1) return int_1;
  127.         /* Speed-up for e.g. 1E-100000 */
  128.  
  129.     if (IsSmallInt(u1) && IsSmallInt(v1))
  130.         return (integer)
  131.             MkSmallInt(dig_gcd(SmallIntVal(u1), SmallIntVal(v1)));
  132.  
  133.     FreezeSmallInt(u1, uu1);
  134.     FreezeSmallInt(v1, vv1);
  135.  
  136.     /* Multi-precision binary gcd algorithm */
  137.  
  138.     {    long k = 0;
  139.         integer t, u, v;
  140.  
  141.         u = (integer) Copy((value) u1);
  142.         v = (integer) Copy((value) v1);
  143.  
  144.  /*B1*/        while (even(u) && even(v)) {
  145.             u = int_half(u);
  146.             v = int_half(v);
  147.             if (++k < 0) {
  148.                 /*It's a number we can't cope with,
  149.                   with too many common factors 2.
  150.                   Though the user can't help it,
  151.                   the least we can do is to allow
  152.                   continuation of the session.
  153.                 */
  154.                 error(MESS(902, "exceptionally large rational number"));
  155.                 k = 0;
  156.             }
  157.         }
  158.  
  159.  /*B2*/        if (even(u)) t = (integer) Copy(u);
  160.         else {
  161.             t = int_neg(v);
  162.             goto B4;
  163.         }
  164.  
  165.         do {
  166.  /*B3*/            do {
  167.                 t = int_half(t);
  168.  B4:                ;
  169.             } while (even(t));
  170.  
  171.  /*B5*/            if (Msd(t) >= 0) {
  172.                 release((value) u);
  173.                 u = t;
  174.             }
  175.             else {
  176.                 release((value) v);
  177.                 v = int_neg(t);
  178.                 release((value) t);
  179.             }
  180.  
  181.  /*B6*/            t = int_diff(u, v);
  182.             /* t cannot be int_1 since both u and v are odd! */
  183.         } while (t != int_0);
  184.  
  185.         release((value) t);
  186.         release((value) v);
  187.  
  188.         while (--k >= 0) u = twice(u);
  189.         
  190.         return u;
  191.     }
  192. }
  193.